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Abstract. This paper settles the computational complexity of 
the problem of integrating a polynomial function / over a rational 
simplex. We prove that the problem is NP-hard for arbitrary poly- 
nomials via a generalization of a theorem of Motzkin and Straus. 
On the other hand, if the polynomial depends only on a fixed num- 
ber of variables, while its degree and the dimension of the simplex 
are allowed to vary, we prove that integration can be done in poly- 
nomial time. As a consequence, for polynomials of fixed total de- 
gree, there is a polynomial time algorithm as well. We conclude 
the article with extensions to other polytopes and discussion of 
other available methods. 



1. Introduction 

Let A be a d- dimensional rational simplex inside M" and let / G 
Q[a;i, . . . , x„] be a polynomial with rational coefficients. We consider 
the problem of how to efficiently compute the exact value of the in- 
tegral of the polynomial / over A, which we denote by f dm. We 
use here the integral Lebesgue measure dm on the affine hull (A) of 
the simplex A, defined below in section 2.1. This normalization of the 
measure occurs naturally in Euler-Maclaurin formulas for a polytope 
P, which relate sums over the lattice points of P with certain inte- 
grals over the various faces of P. For this measure, the volume of the 
simplex and every integral of a polynomial function with rational co- 
efficients are rational numbers. Thus the result has a representation 
in the usual (Turing) model of computation. This is in contrast to 
other normalizations, such as the induced Euclidean measure, where 
irrational numbers appear. 

The main goals of this article are to discuss the computational com- 
plexity of the problem and to provide methods to do the computation 
that are both theoretically efficient and have reasonable performance 
in concrete examples. 

Computation of integrals of polynomials over polytopes is funda- 
mental for many applications. We already mentioned summation over 
lattice points of a polytope. They also make an appearance in recent 
results in optimization problems connected to moment matrices [23]. 
These integrals are also commonly computed in finite element methods, 
where the domain is decomposed into cells (typically simplices) via a 
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mesh and complicated functions are approximated by polynomials (see 
for instance [31]). When studying a random univariate polynomial p{x) 
whose coefficients are independent random variables in certain inter- 
vals, the probability distribution for the number of real zeros of p{x) is 
given as an integral over a polytope [9]. Integrals over polytopes also 
play a very important role in statistics, see, for instance, [25]. Remark 
that among all polytopes, simplices are the fundamental case to con- 
sider for integration since any convex polytope can be triangulated into 
finitely many simplices. 

Regarding the computational complexity of our problem, one can ask 
what happens with integration over arbitrary polytopes. It is very ed- 
ucational to look first at the case when / is the constant polynomial 1, 
and the answer is simply a volume. It has been proved that already 
computing the volume of polytopes of varying dimension is T^P-hard 
[16, 12, 20, 24], and that even approximating the volume is hard [17]. 
More recently in [28] it was proved that computing the centroid of a 
polytope is #P-hard. In contrast, for a simplex, the volume is given by 
a determinant, which can be computed in polynomial time. One of the 
key contributions of this paper is to settle the computational complex- 
ity of integrating a non-constant polynomial over a simplex. Before 
we can state our results let us understand better the input and output 
of our computations. Our output will always be the rational number 
f dm in the usual binary encoding. The d-dimensional input sim- 
plex will be represented by its vertices Si, . . . , s^+i (a ^^-representation) 
but note that, in the case of a simplex, one can go from its represen- 
tation as a system of linear inequalities (an //-representation) to a 
^-representation in polynomial time, simply by computing the inverse 
of a matrix. 

Thus the encoding size of A is given by the number of vertices, the di- 
mension, and the largest binary encoding size of the coordinates among 
vertices. Computations with polynomials also require that one specifies 
concrete data structures for reading the input polynomial and to carry 
on the calculations. There are several possible choices. One common 
representation of a polynomial is as a sum of monomial terms with 
rational coefficients. Some authors assume the representation is dense 
(polynomials are given by a list of the coefficients of all monomials up 
to a given total degree r), while other authors assume it is sparse (poly- 
nomials are specified by a list of exponent vectors of monomials with 
non-zero coefficients, together with their coefficients). Another popu- 
lar representation is by straight-line programs. A straight-line program 
which encodes a polynomial is, roughly speaking, a program without 
branches which enables us to evaluate it at any given point (see [14, 26] 
and references therein). As we explain in Section 2, general straight- 
line programs are too compact for our purposes, so instead we restrict 
to a subclass we call single-intermediate-use (division- free) straight-line 
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programs or SIU straight-line programs for short. The precise defini- 
tion and explanation will appear in Section 2 but for now the reader 
should think that polynomials are represented as fully parenthesized 
arithmetic expressions involving binary operators + and x . 

Now we are ready to state our first result. 

Theorem 1 (Integrating general polynomials over a simplex is hard). 
The following problem is NP-hard. Input: 

(11) numbers d,n ^ N in unary encoding, 

(12) affinely independent rational vectors Si, . . . , s^+i G in binary 
encoding, 

(13) an SIU straight-line program $ encoding a polynomial f G Q[xi, . . . 
with rational coefficients. 

Output, in binary encoding: 

(Oi) the rational number f dm, where A C zs the simplex with 
vertices Si, . . . , s^+i and dm is the integral Lebesgue measure 
of the rational affine subspace (A) . 

But we can also prove the following positive results. 

Theorem 2 (Efficient integration of polynomials of fixed effective 
number of variables). For every fixed number D G N, there exists a 
polynomial-time algorithm for the following problem. 
Input: 

(11) numbers d,n, M E N in unary encoding, 

(12) affinely independent rational vectors Si, . . . , s^^+i G Q" in binary 
encoding, 

(13) a polynomial f G Q[Xi, . . . ,Xd] represented by either an SIU 
straight-line program $ of formal degree at most M, or a sparse 
or dense monomial representation of total degree at most M, 

(14) a rational matrix L with D rows and n columns in binary en- 
coding, the rows of which define D linear forms x 1— >■ on 
M". 

Output, in binary encoding: 

(Oi) the rational number /((^i, x), . . . , (^/j, x)) dm, where A C 
is the simplex with vertices Si, . . . , s^+i and dm is the inte- 
gral Lebesgue measure of the rational affine subspace (A) . 

In particular, the computation of the integral of a power of one linear 
form can be done by a polynomial time algorithm. This becomes false 
already if one considers powers of a quadratic form instead of powers 
of a linear form. Actually, we prove Theorem 1 by looking at powers 
of the Motzkin-Straus quadratic form of a graph. 

Our method relies on properties of integrals of exponentials of linear 
forms. A. Barvinok had previously investigated these integrals and 
their computational complexity (see [3], [5]). 
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As we will see later, when its degree is fixed, a polynomial has a 
polynomial size representation in either the SIU straight-line program 
encoding or the sparse or dense monomial representation and one can 
switch between the three representations efficiently. The notion of for- 
mal degree of an SIU straight-line program will be defined in Section 2. 

Corollary 3 (Efficient integration of polynomials of fixed degree). For 

every fixed number M E N, there exists a polynomial-time algorithm 
for the following problem. Input: 

(11) numbers d,n E N in unary encoding, 

(12) affinely independent rational vectors Si, . . . , s^^+i G Q" in binary 
encoding, 

(13) a polynomial / G Q[ ] represented by either an SIU 
straight-line program $ of formal degree at most M, or a sparse 
or dense monomial representation of total degree at most M. 

Output, in binary encoding: 

(Oi) the rational number f{x.) dm, where A C M" is the sim- 
plex with vertices Si, . . . , s^+i and dm is the integral Lebesgue 
measure of the rational affine subspace (A). 

Actually, we give two interesting methods that prove Corollary 3. 
First, we simply observe that a monomial with total degree M involves 
at most M variables. The other method is related to the polynomial 
Waring problem: we decompose a homogeneous polynomial of total 
degree M into a sum of M-th powers of linear forms. 

In [22] Lasserre and Avrachenkov compute the integral f{x) dm 
when / is a homogeneous polynomial, in terms of the corresponding 
polarized symmetric multilinear form (Proposition 18). We show that 
their formula also leads to a proof of Corollary 3. Furthermore, sev- 
eral other methods can be used for integration of polynomials of fixed 
degree. We discuss them in Section 4. 

This paper is organized as follows: After some preparation in Sec- 
tion 2, the main theorems are proved in Section 3. In Section 4, we 
discuss extensions to other convex polytopes and give a survey of the 
complexity of other algorithms. Finally, in Section 5, we describe the 
implementation of the two methods of Section 3, and we report on a 
few computational experiments. 

2. Preliminaries 
In this section we prepare for the proofs of the main results. 

2.1. Integral Lebesgue measure on a rational afRne subspace of 

M". On itself we consider the standard Lebesgue measure, which 
gives volume 1 to the fundamental domain of the lattice Z". Let L 
be a rational linear subspace of dimension d < n. We normalize the 
Lebesgue measure on L, so that the volume of the fundamental domain 
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Table 1. The representation of {xf + ■ ■ ■ + x"^)^ as a 
straight-line program 



Intermediate 



Comment 



13 
17 







q2 ■ q2 

qi + qs 

X2 



Thus q4 




gs ■ q5 

qi + qe 



Now qj 



X2 I ^2 
1 "I" ' 



qSn-l — 

q3n = qSn-l • qSn-l 

qSn+l = qZn-2 + gSn NoW ^3^+1 = H h 

Q'3n+2 = 1 

•ySn+S = q-in+2 ' ?3n+l 
Q'3n+4 = q^n+Z ' Q'Sn+l 

g3„+fc+2 = gsn+fc+i ■ gsn+i Final result. 



of the intersected lattice L n is 1. Then for any affine subspace 
L + a parallel to L, we define the integral Lebesgue measure dm by 
translation. For example, the diagonal of the unit square has length 1 
instead of ^/2. 

2.2. Encoding polynomials for integration. We now explain our 
encoding of polynomials as SIU straight-line programs and justify our 
use of this encoding. We say that a polynomial / is represented as a 
(division-free) straight-line program $ if there is a finite sequence of 
polynomial functions of Q[a;i, . . . , Xn], namely gi, . . . , g^, the so-called 
intermediate results, such that each qi is either a variable xi, . . . , Xn, an 
element of Q, or either the sum or the product of two preceding poly- 
nomials in the sequence and such that qk = f ■ A straight-line program 
allows us to describe in polynomial space polynomials which otherwise 
would need to be described with exponentially many monomial terms. 
For example, think of of the representation of {x\ + ■ ■ ■ + x'^)^ as mono- 
mials versus its description with only ?,n-\-k-\-2 intermediate results; see 
Table 2. The number of intermediate results of a straight line program 
is called its length. To keep track of constants we define the size of an 
intermediate result as one, unless the intermediate result is a constant 
in which case its size is the binary encoding size of the rational num- 
ber. The size of a straight-line program is the sum of the sizes of the 
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intermediate results. The formal degree of an intermediate result qi is 
defined recursively in the obvious way, namely as if is a constant 
of Q, as 1 if qi is a variable the maximum of the formal degrees 

of the summands if q^ is a sum, and the sum of the formal degrees 
of the factors if is a product. The formal degree of the straight- 
line program $ is the formal degree of the final result qk- Clearly the 
total degree of a polynomial is bounded by the formal degree of any 
straight-line program which represents it. 

A favorite example to illustrate the benefits of a straight-line pro- 
gram encoding is that of the symbolic determinant of an n x n matrix. 
Its dense representation as monomials has size 0(n!) but it can be com- 
puted in 0{n^) operations by Gaussian elimination. See the book [14] 
as a reference for this concept. 

From a monomial representation of a polynomial of degree M and 
n variables it is easy to encode it as a straight-line program: first, by 
going in increasing degree we can write a straight-line program that 
generates all monomials of degree at most M in n variables. Then for 
each of them compute the product of the monomial with its coefficient 
so the length doubles. Finally successively add each term. This gives 
a final length bounded above by four times the number of monomials 
of degree at most M in n variables. 

Straight-line programs are quite natural in the context of integra- 
tion. One would certainly not expand {x\ + ■ ■ ■ + to carry on 
numeric integration when we can easily evaluate it as a function. More 
importantly, straight-line programs are suitable as an input and output 
encoding and data structure in certain symbolic algorithms for com- 
putations with polynomials, like factoring; see [14]. Since straight-line 
programs can be very compact, the algorithms can handle polynomials 
whose input and output encodings have an exponential size in a sparse 
monomial representation. 

However, a problem with straight-line programs is that this input 
encoding can be so compact that the output of many computational 
questions cannot be written down efficiently in the usual binary en- 
coding. For example, while one can encode the polynomial x"^ with 
a straight-line program with only k + 1 intermediate results (see Ta- 
ble 2), when we compute the value of x"^ for x = 2, or the integral 
Jq x'^'' dx = 2^''"*"^/ (2'^ + 1), the binary encoding of the output has a 
size of 9(2^). Thus the output, given in binary, turns out to be ex- 
ponentially bigger than the input encoding. We remark that the same 
difficulty arises if we choose a sparse input encoding of the polynomial 
where not only the coefficients but also the exponent vectors are en- 
coded in binary (rather than the usual unary encoding for the exponent 
vectors). 
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Table 2. The representation of a straight-line program 
for , using iterated squaring 



Intermediate Comment 





= X 






= Qi ■ 




^3 


= q2 ' 






= Qk 


■ qk Final result. 



Table 3. The representation of a single-intermediate- 
use straight-line program for x"^ ; note that the iterated 



squaring method cannot be 


used 


Intermediate 


Comment 


qi= X 




q2 = X 




q3 = qi- q2 


Now ^3 = x^, and qi and q2 




cannot be used anymore. 


qA = X 




q5 = qs- q4 


Thus gs = x^. 


q2k+l_2 = X 




q2k+i_i = q2k+i-s ■ q2h+i-2 


Final result. 



This motivates the following variation of the notion of straight- 
line program: We say a (division-free) straight-line program is single- 
intermediate-use, or SIU for short, if every intermediate result is used 
only once in the definition of other intermediate results. (However, 
the variables be used arbitrarily often in the definition 

of intermediate results.) With this definition, all ways to encode the 
polynomial require at least 2*^ multiplications. An example SIU 
straight-line program is shown in Table 3. Clearly single-intermediate- 
use straight-line programs are equivalent, in terms of expressiveness 
and encoding complexity, to fully parenthesized arithmetic expressions 
using binary operators + and x . 

2.3. Efficient computation of truncated product of an arbi- 
trary number of polynomials in a fixed number of variables. 

The following result will be used in several situations. 

Lemma 4. For every fixed number D G N, there exists a polynomial 
time algorithm for the following problem. 



8 V. BALDONI, N. BERLINE, J. A. DE LOERA, M. KOPPE, AND M. VERGNE 



Input: a number M in unary encoding, a sequence of k polynomials 
Pj e Q[Xi, . . . , Xu] of total degree at most M, in dense monomial 
representation. 

Output: the product Pi ■ ■ ■ truncated at degree M. 

Proof. We start with the product of the first two polynomials. We 
compute the monomials of degree at most M in this product. This takes 
0(M^^) elementary rational operations, and the maximum encoding 
length of any coefficient in the product is also polynomial in the input 
data length. Then we multiply this truncated product with the next 
polynomial, truncating at degree M, and so on. The total computation 
takes 0(/cM^^) elementary rational operations. □ 

3. Proofs of the main results 

Our aim is to perform an efficient computation of f dm where 
A is a simplex and / a polynomial. We will prove first that this is 
not possible for / of varying degree under the assumption that P 7^ 
NP. More precisely, we prove that, under this assumption, an efficient 
computation of dm is not possible, where Q is a quadratic form 

and M is allowed to vary. 

In the next subsection we present an algorithm to efficiently compute 
the integral f dm in some particular situations, most notably the 
case of arbitrary powers of linear forms. 

3.1. Hardness for polynomials of non-fixed degree. For the proof 
of Theorem 1 we need to extend the following well-known result of 
Motzkin and Straus [27]. In this section, we denote by A the (n — 1)- 
dimensional canonical simplex {x e M" : Xi > 0, Yl^=i^i — 
we denote by dm the Lebesgue measure on the hyperplane { x G : 
^27=1 = 1}, normalized so that A has volume 1. For a function / on 
A, denote as usual ||/||^ = max^gA |/(x)| and \\f\\^ = (/^ \ff dm)^/f, 
for p > 1. Recall that the clique number of a graph G is the largest 
number of vertices of a complete subgraph of G. 

Theorem 5 (Motzkin-Straus). Let G be a graph with n vertices and 
clique number uj{G) . Let Qg{'^) be the Motzkin-Straus quadratic form 
lE(M)ei?(G)^*^r Then ||Qg(x)|L = " 

Our first result might be of independent interest as it shows that in- 
tegrals of polynomials over simplices can carry very interesting combi- 
natorial information. This result builds on the theorem of Motzkin and 
Straus, using the proof of the well-known relation \\f\\^ = limp_,oo 11/ Hp- 
Lemma 6. Let G be a graph with n vertices and clique number uj{G). 
Let (5g(x) be the Motzkin-Straus quadratic form. Then forp > 4(e — 
l)n^ ln(32n^), the clique number u{G) is equal to [ i_2||q 11 ] ■ 
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To prove Lemma 6 we will first prove the following intermediate 
result. 

Lemma 7. For e > we have 

i\\QG\L-e){'-r-'^/'<\\Qal<\\Qa\\^. 

Proof. The right-hand side inequality follows from the normalization 
of the measure, as |(5(x)| < ||Qg|Ioo' x G A. 

In order to obtain the other inequality, we use Holder's inequality 
/aI/A'I dm < ll/llp llS'llg, where q is such that i + i = 1. For any 
(say) continuous function / on A, let us denote by A{f,e) the set 
{xG A : |/(x)| > ll/lloo"^}' ^iid take for g the characteristic function 
of A(/,e). We obtain 

,-e)(volA(/,e))Vp<||/|| . (1) 



Let a be a point of A where the maximum of Qg is attained. Since 
^ = J2ii,j)eE{G) W6 know that < ^ < 1 for x G A. Since A is 
convex, we conclude that for any x G A, 



< gG(a) - Qg{^) < 5^ - 



1=1 



Thus A{Qg,£) contains the set = { x G A : XliLi 1*^* ~ < 
We claim that vo^C^) > (f)'"-^ . This claim proves the left inequality 
of the lemma when we apply it to (1). 

Consider the dilated simplex -^^^A and the translated set = 

+ i+e/2 ^- Clearly P^ is contained in A. Moreover, for x G P^, 
we have X]"=i ^ ^ i+e/2 — ^' hence Pg is contained in Ce- Since 
vol(A) = 1 for the normalized measure, the volume of Pg is equal to 
(t+^)"~^- Hence vo^P^) > This finishes the proof. □ 

Proof of Lemma 6. In the inequalities of Lemma 7, we substitute the 
relation ||Qg|Ioo ~ ~ Ij(by)' S^'^^^ ^y Motzkin-Straus's theorem 
(Theorem 5). We obtain 

(^(i-^)-^)(^/4)"^<IIQgII <i(i ^ ^ 



Let us rewrite these inequalities as 

T^^^^-'-G)< . (2) 

We only need to prove that for e = and p > 4(e — l)n^ ln(32n^) we 
have 

< m := — ^ ] < 1. (3) 

(£/4) — 
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Let us write 

^ = WQgI {—Kei - 1) = WQgI mn')"^ 

Thus L{p) in Equation 3 becomes now 

1 

'HQg 

Since HQcUp < |(1 — Zj(g)^ — h have a bound for 6p 
0<Sp< ^ii32n^)'^ - 1). 

Let A= = (32n2)"-i. Since we assumed p > 4(e-l)n3 ln(32ra2) 

we have o'< ^ < 1, hence < A^/p - 1 < (e - 1)^^^^. We obtain 

o<. <£^("-i)i°gP2"-)< 1 



^ - 2 p - 8n2 ■ 

Since uj{G) < n, we have 1 — 2 HQdlp > 1/n. Hence we have 

2{5, + s) 



1-2\\Qg\\p - 2n- 2- 

Finally for any number < a < 1/2 we have y~ < 1 + 2a, hence 
applying this fact to Equation 4 with a = we get 

This proves Equation 3 and the lemma. □ 

Proof of Theorem 1. The problem of deciding whether the clique num- 
ber uj{G) of a graph G is greater than a given number i^' is a well-known 
NP-complete problem [18]. From Lemma 6 we see that checking this 
is the same as checking that for p = 4(e — l)n^ ln(32n^) the integral 
part of J^^QgY dm is less than K^. Note that the polynomial Qd'x.y 
is a power of a quadratic form and can be encoded as a SIU straight- 
line program of length 0(n^logn ■ \E{G)\). If the computation of the 
integral f dm of a polynomial / could be done in polynomial time 
in the input size of /, we could then verify the desired inequality in 
polynomial time as well. □ 

3.2. An extension of a formula of Brion. In this section, we obtain 
several expressions for the integrals dm and ■ ■ ■ i'j^^ dm, 

where A C M" is a simplex and i, ij are linear forms on M". The 
first formula, (5) in Lemma 8, is obtained by elementary iterated in- 
tegration on the standard simplex. It leads to a computation of the 
integral ■ ■ ■ i'^'^ dm in terms of the Taylor expansion of a certain 



HOW TO INTEGRATE A POLYNOMIAL OVER A SIMPLEX 



11 



analytic function associated to A (Corollary 11), hence to a proof of 
the complexity result of Theorem 2. 

In the case of one linear form £ which is regular, we recover in this 
way the "short formula" of Brion as Corollary 12. This result was first 
obtained by Brion as a particular case of his theorem on polyhedra [13]. 

Lemma 8. Let A be the simplex that is the convex hull of{d+l) affinely 
independent vertices Si, S2, . . . , s^+i in M", and let I he an arbitrary 
linear form on M". Then 

/ e dm = rf! vol(A, dm) (Ikl+dV. ' 

where |k| = J2'j=i kj- 

Proof. Using an affine change of variables, it is enough to prove (5) 
when A is the d-dimensional standard simplex A^t C M"^ defined by 

A,i=|xGM^:Xi>0,^a;i<l|. 

The volume of A^t is equal to ^. In the case of A^f, the vertex is 
the basis vector e^- for 1 < j < and s^+i = 0. Let (£, x) = J2'j=i ^j^j- 
Then (5) becomes 

fci fcrf 



/ gQixiH ha^xa _ 



. (|k|+d)! 

We prove it by induction on d. For = 1, we have 

e-- dx = - - = V 

k>0 ^ ' 







Let d> 1. We write 

gOixiH l-a^Xi 



Xl-\ \-Xci-l<l~Xd 

By the induction hypothesis and an obvious change of variables, the 
inner integral is equal to 

(1-x.)^"^ (1-^.^"^' ' 



(Ikl 1)!' 



The result now follows from the relation 

dx = ^ 



^ il-xY _ , a" 



_^ax. 



p\ ^(A;+p+l)!' 
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□ 

Remark 9. Let us replace I by ti in (5) and expand in powers of t. 
We obtain the following formula. 

I t'&m = vol(A, dm) J] {i,s,)'^ ■ ■ ■ {i,s,^,)''^^K 

(6) 

This relation is a particular case of a result of Lasserre and Avrachenkov, 
Proposition 18, as we will explain in section 4.3 below. 

Theorem 10. Let A be the simplex that is the convex hull of [d + 1) 
affinely independent vertices Si, S2, . . . , s^+i in M". 



y ^M iM + d)\ I ^M^^^ ^,^^1^^ dm)— -^J^ 



-• (7) 

)) 



Proof. We apply Formula (6). Summing up from M = to cxd, we 
recognize the expansion of the right-hand side of (7) into a product of 
geometric series: 



^ Ml L 



d\ vol(A, dm) ^ J2 ^i)'' ■ ■ ■ S'^+i)''^'- 

Mm kGN''+i|,k|=M 

□ 

Theorem 10 has an extension to the integration of a product of pow- 
ers of several linear forms. The following formula is implemented in 
our Maple program duality .mpl, see Table 8 

Corollary 11. Let ii, . . . ,£0 be D linear forms on M". We have the 
following Taylor expansion: 

2^/1 ■■■'o rf,vol(A,dm)AM,!...M^!^"^- 



n-=i(i-ti(^i,s.) tn{iD,s,)y 

Proof. Replace t£ with tiii + ■ — h t£)i£) in (7) and take the expansion 
in powers □ 

From Theorem 10, we obtain easily the "short formula" of Brion, in 
the case of a simplex. 

Corollary 12 (Brion). Let A be as in the previous theorem. Let i be a 
linear form which is regular w.r.t. A, i.e., {i,Si) 7^ Sj) for any pair 
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i ^ j ■ Then we have the following relations. 

/ i^' dm = d\ vol( A, dm) , J ( V .^^^ A . (9) 



(10) 



Proof. We consider the right-hand side of (7) as a rational function of 
t. The poles t = Sj) are simple precisely when £ is regular. In this 
case, we obtain (9) by taking the expansion into partial fractions. The 
second relation follows immediately by expanding e^. □ 

When i is regular, Brion's formula is very short, it is a sum of c? + 1 
terms. When £ is not regular, the expansion of (7) into partial fractions 
leads to an expression of the integral as a sum of residues. Let K C 
{1, . . . , + 1} be an index set of the different poles t = 1/ (£, Sfc), and 
for G -ft' let mfc denote the order of the pole, i.e., 

mfc = #{ i G {1, . . . , + 1} : (£, s,) = {£, Sk) }• 

With this notation, we have the following formula, which is imple- 
mented in our Maple program waring. mpl, see Tables 5 and 6. 

Corollary 13. 



(111 



Remark 14. It is worth remarking that Corollaries 12 and 13 can 
be seen as a particular case of the localization theorem in equivariant 
CO homology (see for instance [8]), although we did not use this fact and 
instead gave a simple direct calculation. In our situation, the variety is 
the complex projective space CP"^, with action of a rf-dimensional torus, 
such that the image of the moment map is the simplex A. Brion's 
formula corresponds to the generic case of a one-parameter subgroup 
acting with isolated fixed points. In the degenerate case when the 
set of fixed points has components of positive dimension, the polar 
parts in (11) coincide with the contributions of the components to the 
localization formula. 

A formula equivalent to Corollary 13 appears already in [3] (3.2). 



3.3. Polynomial time algorithm for polynomial functions of a 
fixed number of linear forms. 
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Proof of Theorem 2. We now present an algorithm which, given a poly- 
nomial of the particular form x), ■ ■ ■ , {£d,x.)) where / is a poly- 
nomial depending on a fixed number D of variables, and = 
LjiXi -!-■■■ + LjnXn, for j = 1, . . . , D, are linear forms on M", computes 
its integral on a simplex, in time polynomial on the input data. This 
algorithm relies on Corollary 11. 

The number of monomials of degree M in D variables is equal to 
i^^D-i^)' Therefore, when D is fixed, the number of monomials of 
degree at most M in D variables is O(M^). When the number of 
variables D of a straight-line program $ is fixed, it is possible to com- 
pute a sparse or dense representation of the polynomial represented 
by $ in polynomial time, by a straight-forward execution of the pro- 
gram. Indeed, all intermediate results can be stored as sparse or 
dense polynomials with O(M^) monomials. Since the program $ is 
single-intermediate-use, the binary encoding size of all coefficients of 
the monomials can be bounded polynomially by the input encoding 
size. Thus it is enough to compute the integral of a monomial, 

/ (^i,x)*^^---(£z),x)*^-dm. (12) 

J A 

From Corollary 11, it follows that 



(|M| +d)! 



' '"^^ dm 



rf!vol(A, dm) 7a Mi\---Md\ 
is the coefficient of t*^^ ■ ■ ■ in the Taylor expansion of 



n?=l(l-^l(^l>S.) tD(^D,S,,))' 

Since D is fixed, this coefficient can be computed in time polynomial 
with respect to M and the input data, by multiplying series truncated 
at degree |M|, as explained in Lemma 4. 

Finally, vol(A, dm) needs to be computed. If A = convjsi, . . . , s^+i} 
is full- dimensional {d = n), we can do so by computing the determinant 
of the matrix formed by difference vectors of the vertices: 

vol(A, dm) = ^ |det(si - s„+i, . . . , s„ - s„+i)| . 

If A is lower- dimensional, we first compute a basis B G Z"^^"^ of the 
intersection lattice A = lin(A) (1 Z". This can be done in polynomial 
time by applying an efficient algorithm for computing the Hermite nor- 
mal form [19]. Then we express each difference vector Vj = Sj — 8,^+1 G 
lin(A) for z = 1, . . . , d using the basis B as Vj = BV^, where G Q"'. 
We obtain 

vol(A, dm) = ^|det(v;,...,v:,)|. 
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thus the volume computation is reduced to the calculation of a deter- 
minant. This finishes the proof of Theorem 2. □ 

3.4. Polynomial time algorithms for polynomials of fixed de- 
gree. In the present section, we assume that the total degree of the 
input polynomial / we wish to integrate is a constant M. 

Proof of Corollary 3. First of all, when the formal degree M of a straight- 
line program $ is fixed, it is possible to compute a sparse or dense 
representation of the polynomial represented by $ in polynomial time, 
by a straight-forward execution of the program. Indeed, all intermedi- 
ate results can be stored as sparse or dense polynomials with 0{n'^^) 
monomials. Since the program $ is single-intermediate-use, the bi- 
nary encoding size of all coefficients of the monomials can be bounded 
polynomially by the input encoding size. 

Now, the key observation is that a monomial of degree at most M 
depends effectively on D < M variables Xj^, . . . , Xj^, thus it is of the 
form 

'^1 ' ' ' '^D -I 

where the linear forms •^j(x) = Xi. are the coordinates that effec- 
tively appear in the monomial. Thus, Corollary 3 follows immediately 
from Theorem 2. This method is implemented in our Maple program 
duality. mpl, see Tables 8 and 12. □ 

Remark 15. The relations in Corollary 12 can be interpreted as equal- 
ities between meromorphic functions of i. The right-hand side is a sum 
of meromorphic functions whose poles cancel out, so that the sum is ac- 
tually analytic. We derive from this another polynomial time algorithm 
for computating the integral 

/ x'^^-'-x^' dm. 

More precisely, let us write (£, x) = y\X\ + ■ ■ ■ + ydXd- Then the integral 
x™^ ■ ■ ■ x™'* dm is the coefficient of ^^^'''^^i in the Taylor expansion 
of gyi^i+ - !/d3^d ^j^^ compute it by taking the expansion of each 
of the terms of the right hand-side of Equation (10) into an iterated 
Laurent series with respect to the variables Xi, . . . , x^. This method is 
implemented in our Maple program iterated-laurent .mpl, see Ta- 
bles 11 and 14. 

In the following, we give another proof of Corollary 3, based on 
decompositions of polynomials as sums of powers of linear forms. 

Alternative 'proof of Corollary 3. From Corollaries 12 and 13, we derive 
another efficient algorithm, as follows. The key idea now is that one 
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can decompose the polynomial / as a sum / := Ylii^f-^f with at most 
2^ terms in the sum. We use the well-known identity 



IMI 



|M|! 



(13) 

where |M| = Mi H h M„ < M. 

In the implementation of this method, we may group together pro- 
portional linear forms. The number F{n, M) of primitive vectors (pi, . . . , p„) 
which appear in the decomposition of a polynomial of total degree < M 
is given by the following closed formula. ^. 

Lemma 16. Let 

F{n,M) = 

Card({(pi, . . . ,p.) e N",gcd(pi, . . . ,p„) = 1, 1 < < M}). 

i 

Then 

F{n,M) = J2Mi[ (14) 
d=i ^ ^ 
where fi{d) is the Mobius function. 
When M is fixed and n ^ oo we have 

M 

i^(^,M) = ^ + 0(n^-i). 

Proof. Let G(n, M) = Card({(pi, . . . ,p„) G N", 1 < Y.Pi < M}). By 
grouping together the vectors {pi, ■ ■ ■ ,Pn) with a given gcd d, we obtain 

^ M 
G{n,M)=J2F{n, [-j])- 

d=l 

Moreover the number of all integral vectors [pi, . . . ,pn) € N" such that 
Yli Pi ^ M is equal to the binomial coefficient ("'^^) • When we omit 
the zero- vector we obtain 

Then we obtain (14) by applying the second form of Mobius inversion 
formula. The asymptotics follow easily from (14). □ 

Thus Formula (13), together with Corollary 13, give another poly- 
nomial time algorithm for integrating a polynomial of fixed degree. It 
is implemented in our Maple program waring. mpl, see Tables 6 and 
10. □ 



"'^Lemma 16 was kindly supplied by Christophe Margerin. 
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The problem of finding a decomposition with the smallest possible 
number of summands is known as the polynomial Waring problem. 
Alexander and Hirschowitz solved the generic problem (see [1], and 
[11] for an extensive survey). 

Theorem 17. The smallest integer r{M,n) such that a generic homo- 
geneous polynomial of degree M inn variables is expressible as the sum 
ofr{M,n) M-th powers of linear forms is given by 



with the exception of the cases r(3,5) = 8, r(4, 3) = 6, r(4, 4) = 10, 
r(4, 5) = 15, and M = 2, where r(2, n) = n. 

An algorithm for decomposing a given polynomial into the smallest 
possible number of powers of linear forms can be found in [10]. 

In the extreme case, when the polynomial / happens to be the power 
of one linear form i, one should certainly avoid applying the above 
decomposition formula to each of the monomials of /. We remark that, 
when the degree is fixed, we can decide in polynomial time whether a 
polynomial /, given in sparse or dense monomial representation, is a 
power of a linear form i, and, if so, construct such a linear form. 



4. Other algorithms for integration and extensions to 



We conclude with a discussion of how to extend integration to other 
polytopes and a review of the complexity of other methods to integrate 
polynomials over polytopes. 

4.1. A formula of Lasserre— Avrachenkov. Another nice formula 
is the Lasserre-Avrachenkov formula for the integration of a homoge- 
neous polynomial [22] on a simplex. As we explain below, this yields a 
polynomial-time algorithm for the problem of integrating a polynomial 
of fixed degree over a polytope in varying dimension, thus providing an 
alternative proof of Corollary 3. 

Proposition 18 ([22]). Let H be a symmetric multilinear form defined 
on (M.'^)'^^ . Let Si, S2, . . . , s^+i be the vertices of a d- dimensional simplex 
A. Then one has 



r{M,n) 



- (n+M~l\ - 
\ M ) 



n 



other polytopes 
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Remark 19. By reindexing the summation in (15), as H is symmetric, 
we obtain 



/ i?(x, X, . . . , x)(ix 



vol (A) 

(M+d\ 
\ M ) ki 



^ i7(Si, . . . ,Si, . . . ,Sd+i, . . . ,Sd+i), (16) 

where Si is repeated ki times, S2 is repeated k2 times, etc. When H 
is of the form H = nf=i(^'^«)' ^ single linear form i, then (16) 
coincides with Formula (6) in Remark 9. 

Now any polynomial / which is homogeneous of degree M can be 
written as /(x) = Hf{:K, x, . . . , x) for a unique multilinear form Hf. If 

/ = then Hf = YlfLii^^^i)- Thus for fixed M the computation of 
Hf can be done by decomposing / into a linear combination of powers 
of linear forms, as we did in the proof of Corollary 3. Alternatively one 
can use the well-known polarization formula, 

1 

i7j(xi,...,XM) = ^ig^ ^ eie2---eMf(j2eiX,^, (17) 
Thus from (15) we get the following corollary. 

Corollary 20. Let f be a homogeneous polynomial of degree M in 
d variables, and let Si,S2, . . . ,8^4.1 be the vertices of a d-dimensional 
simplex ^. Then 

fiy) dy = 

vol (A) 



eie2- ■■SMfi^SkSi^. (18) 

l<ji<j2<---<i]\f<d+l ee{±i}M k=l 



We remark that when we fix the degree M of the homogeneous poly- 
nomial /, the length of the polarization formula (thus the length of the 
second sum in (18) is a constant. The length of the first sum in (18) is 
0(n^^). Thus, for fixed degree in varying dimension, we obtain another 
polynomial-time algorithm for integrating over a simplex. 

4.2. Traditional conversion of the integral as iterated univari- 
ate integrals. Let P C M°' be a full-dimensional polytope and / a 
polynomial. The traditional method we teach our calculus students 
to compute multivariate integrals over a bounded region requires them 
to write the integral Jp / dm is a sum of sequences of one-dimensional 
integrals 

^ / / ■ ■ ■ / / dxi, dxi^ . . . dxi^ (19) 

7=1 "'"ij "'C'2j Ja.dj 
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for which we know the hmits of integration ajj, bij exphcitly. The prob- 
lem of finding the hmits of integration and the sum has interesting com- 
plexity related to the well-known Fourier-Motzkin elimination method 
(see Chapter One in [30] for a short introduction). 

Given a system of linear inequalities Ax < b, describing a polytope 
P C M'^, Fourier-Motzkin elimination is an algorithm that computes a 
new larger system of inequalities Ax < h with the property that those 
inequalities that do not contain the variable Xd describe the projection 
of P into the hyperplane Xd = 0- We will not explain the details, but 
Fourier-Motzkin elimination is quite similar to Gaussian elimination 
in the sense that the main operations necessary to eliminate the last 
variable Xd require to rearrange, scale, and add rows of the matrix 
{A,h), but unlike Gaussian elimination, new inequalities are added to 
the system. 

It was first observed by Schechter [29] that Fourier-Motzkin elim- 
ination provides a way to generate the traditional iterated integrals. 
More precisely, let us call Pd the projection of P into the the hyper- 
plane Xd = 0. Clearly when integrating over a polytopal region we 
expect that the limits of integration will be affine functions. From 
the output of Fourier-Motzkin Ax < h, we have x G P if and only if 
(xi, . . . , Xd-i) e Pd and for the first k + r inequalites of the system 

d-l 

Xd <h-^ aijXj = A^"(xi, . . . , Xd-i) 

5=1 

for z = 1, . . . as well as 

d-l 

Xd > h+i - ^ ak+ijXj, = A-(xi, . . . , Xd~i) 
i=i 

for i = 1, . . . ,r. Then, if we define 

m(xi, ...,Xd)= max{A]{xi, . . . , Xd-i), j = 1, . . . , r} 

and 

M{xi, ...,Xd) = min{A^(xi, . . . , Xd-i), j = 1, . . . , r}, 
we can write 

/ /(x) dm = / / /(x) dxi dx2 ■ ■ ■ dxd 

Jp Jp^ Jrn 



Finally the convex polytope Pd can be decomposed into polyhedral 
regions where the functions m, M become simply affine functions from 
among the list. Since the integral is additive we get an expression 

/ /(x) dm = / / /(x) dxi dx2 ■ ■ ■ dxd- 

Jp^ M 
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Finally by repeating the elimination of variables we recover the full 
iterated list in (19). As it was observed in [29], this algorithm is unfor- 
tunately not efficient because the iterated Fourier-Motzkin elimination 
procedure can produce exponentially many inequalities for the descrip- 
tion of the projection (when the dimension d varies). Thus the number 
of summands considered can in fact grow exponentially. 

4.3. Two formulas for integral computation. We would like to 
review two formulas that are nice and could speed up computation in 
particular cases although they do not seem to yield efficient algorithms 
just on their own. 

First, one may reduce the computation of Jpfdm to integrals on 
the facets of P, by applying Stokes formula. We must be careful to use 
a rational Lebesgue measure on each facet. As shown in ([4]), we have 
the following result. 

Theorem 21. Let {-Fi}i=i,...,m be the set of facets of a full- dimensional 
polytope P C M". For each i, let Uj be a rational vector which is 
transverse to the facet Fi and pointing outwards P and let dfii be the 
Lebesgue measure on the affine hull of Fi which is defined by contracting 
the standard volume form of M" with Uj . Then 

Jp (a, y) ^ Jf, 

for all a G C" and y G M" such that (y, a) 7^ 0. 

It is clear that, by considering the expansion of the analytic function 
jp e^^'^"^ dx, we can again obtain an analogous result for polynomials. 
An alternative proof was provided by [21]. The above theorem, how- 
ever, does not necessarily reduce the computational burden because, 
depending on the representation of the polytope, the number of facets 
can be large and also the facets themselves can be complicated poly- 
topes. Yet, together with our results we obtain the following corollary 
for two special cases. 

Corollary 22. There is a polynomial-time algorithm for the following 
problem. Input: 

(11) the dimension n E'R in unary encoding, 

(12) a list of rational vectors in binary encoding, namely 

(i) either vectors (hi, Hi q), . . . , (h^, h^fi) G Q""*"-^ that describe 
the facet- defining inequalities (hj,x) < hi Q of a simplicial 
full- dimensional rational polytope P , 

(ii) or vectors Si, . . . , Sjy G Q"' that are the vertices of a simple 
full- dimensional rational polytope P , 

(13) a rational vector a G Q" in binary encoding, 

(14) an exponent M ^ N in unary encoding. 
Output, in binary encoding. 
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(Oi) the rational number 

/(x) dm where /(x) = (a, x)^ 

and where dm is the standard Lebesgue measure on M". 

Proof. In the case (i) of simplicial polytopes P given by facet-defining 
inequalities, we can use linear programming to compute in polynomial 
time a ^^-representation for each simplex Fi that is a facet of P. By 
applying Theorem 21 with ta in place of a and extracting the coefficient 
of t'^ in the Taylor expansion of the analytic function t Jp(ta), we 
obtain the formula 

/p<«- = (M+l')(y.a) /, 

which holds for all y G with (y, a) 7^ 0. It is known that a suit- 
able y 6 Q" can be constructed in polynomial time. The integrals 
on the right-hand side can now be evaluated in polynomial time using 
Theorem 2. 

In the case (ii) of simple polytopes P given by their vertices, we make 
use of the fact that a variant of Brion's formula (10) actually holds for 
arbitrary rational polytopes. For a simple polytope P, it takes the 
following form. 

where N{si) denotes the set of vertices adjacent to Sj in P, and A, = 
|det(sj — Sj)jg7v(s;)|. The right-hand side is a sum of rational functions 
of i, where the denominators cancel out so that the sum is actually 
polynomial. If i is regular, that is to say {i, Sj — Sj) ^ for any i and 
j e N(si), then the integral can be computed by (20) which is a very 
short formula. However it becomes difficult to extend the method which 
we used in the case of a simplex. Instead, we can do a perturbation. In 
(20), we replace i hj i + ei', where i' is such that i + ei' is regular for 
e 7^ 0. The algorithm for choosing is bounded polynomially. Then 
we do expansions in powers of e as explained in Lemma 4. □ 

4.4. Triangulation of arbitrary polytopes. It is well-known that 
any convex polytope can be triangulated into finitely many simplices. 
Thus we can use our result to extend the integration of polynomials 
over any convex polytope. The complexity of doing it this way will 
directly depend on the number of simplices in a triangulation. This 
raises the issue of finding the smallest triangulation possible of a given 
polytope. Unfortunately this problem was proved to be NP-hard even 
for fixed dimension three (see [15]). Thus it is in general not a good idea 
to spend time finding the smallest triangulation possible. A cautionary 
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remark is that one can naively assume that triangulations help for non- 
convex polyhedral regions, while in reality it does not because there 
exist nonconvex polyhedra that are not triangulable unless one adds 
new points. Deciding how many new points are necessary is an NP- 
hard problem [15]. 

5. Implementation and computational experiments 

We have written Maple programs to perform some initial experiments 
with the three methods described in Section 3.4. The programs are 
available at [2].^ 

5.1. Integration of a power of a linear form and decomposition 
of polynomials into powers of linear forms. 

5.1.1. Decomposition of polynomials into powers of linear forms. Ta- 
ble 4 shows the number F{n, M) of primitive linear forms {pi, . . . ,Pn) 
which may appear in the decomposition (13) of a polynomial of total 
degree < M. This number is computed using the closed formula (14). 

Table 4. Decomposition of polynomials into powers of 
linear forms 



n 










Degree M 










1 


2 


5 


10 


20 


30 




40 


50 


2 


2 


3 


11 


33 


129 


279 




491 


775 


3 


3 


6 


40 


205 


1381 


4306 




9880 


18970 


4 


4 


10 


103 


831 


9373 


41373 


122349 


286893 


5 


5 


15 


221 


2681 


49586 


305836 


1.2 


• 10^ 


3.3 • 10*^ 


8 


8 


36 


1226 


42271 


3.1 • 10*^ 


4.8 • 10^ 


3.7 


•10^ 


1.9 • 109 


10 


10 


55 


2917 


181413 


3.0 • 10"^ 


8.4 • 10^ 


1.0 • 


IQio 


7.5 • 10i° 


15 


15 


120 


15338 


3.3 • 10^ 


3.2 • 10^ 


3.4 • 10" 


1.2 • 


1013 


2.1 • IQi'' 


20 


20 


210 


52859 


3.0 • 10^ 


1.4 • lO^i 


4.7 • 1013 


4.2 • 


1015 


1.6 • lOi^ 


30 


30 


465 


324076 


8.5 • 10*^ 


4.7-1013 


1.2 • lOi'^ 


5.5 • 


1019 


8.9 • 10^1 


40 


40 


820 


1.2 • 10^ 


1.0 • 10^° 


4.2 • lOi^"^ 


5.5 • 1019 


1.1 • 


1023 


6.0 • 10^5 


50 


50 


1275 


3.5 • 10^ 


7.5 • 10^0 


1.6 • lOi^ 


8.9 • 10^1 


6.0 • 


1025 


1.0 ■ 1029 



AH algorithms are implemented in the files waring.mpl, 
iteratedJLaurent .mpl, duality .mpl; all tables with random examples are 
created using procedures in exEimples .mpl and tables .mpl. 
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5.1.2. Integration of a power of a linear form over a simplex. We have 
written a Maple program which implements the method of Corollary 13 
for the efficient integration of a power of one linear form over a simplex, 
dm."' In a computational experiment, for a given dimension n and 
degree M, we picked random full-dimensional simplices A and random 
linear forms ^ and used the Maple program to compute the integral. 
Table 5 shows the computation times.* 

Table 5. Integration of powers of linear forms over simplices 



Degree M 



n 


2 


10 


20 


50 


100 


300 


1000 


10 


0.0 


0.0 


0.0 


0.1 


0.0 


0.0 


0.0 


20 


0.1 


0.2 


0.1 


0.2 


0.1 


0.2 


0.2 


50 


1.0 


1.4 


1.4 


1.6 


1.6 


1.6 


1.7 


100 


5.1 


8.4 


8.7 


9.2 


9.5 


10.0 


11.0 


200 


36 


71 


84 


88 


97 


110 


120 


300 


150 


320 


400 


470 


520 


530 




400 


500 














1000 

















5.1.3. Integration of a monomial over a simplex by decomposition as 
sum of powers of linear forms. Next, we tested the algorithm which 
computes the integral of a monomial over a simplex A, by de- 
composing it as a sum of powers of linear forms. This algorithm was 
discussed in Section 3.4. In our experiments, for given dimension n and 
total degree M, we picked 50 combinations of a random simplex A of 
dimension n and a random exponent vector M = (Mi, . . . , M„) with 

Er=iM, = M. 

First we decompose a given monomial into a sum of powers of linear 
forms, then we integrate each summand using the Maple procedure 
discussed above. ^ Table 6 shows the minimum, average, and maximum 
computation times. 



The integration is done by the Maple procedure integraLpowerJinear^form in 
waring.mpl. 

*A11 experiments were done with Maple 12 on Sun Fire V440 machines with 
UltraSPARC-IIIi processors running at 1.6 GHz. The computation times are given 
in CPU seconds. All experiments were subject to a time limit of 600 seconds per 
example. 

^This method is implemented in the Maple procedure integraLvia_waring in 
waring.mpl. 
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Table 6. Integration of a random monomial of pre- 
scribed degree by decomposition into a sum of powers 
of linear forms 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 


100 


200 


300 


2 




















0.1 


0.3 


0.8 


6.2 


0.3 













0.1 


0.2 


0.5 


1.0 


1.5 


8.3 


38 


90 







0.1 


0.1 


0.1 


0.4 


0.8 


1.9 


2.6 


12 


59 


174 


3 

















0.4 


0.5 


0.9 


4.6 

















0.2 


1.1 


3.0 


7.4 


17 


173 











0.2 


0.1 


0.4 


2.2 


7.0 


19 


37 


440 






4 











0.1 


0.6 


1.5 


4.8 


9.7 
















0.1 


0.4 


3.6 


15 


52 


135 













0.2 


0.2 


0.9 


8.7 


44 


149 


404 








5 











0.1 


0.3 


4.5 




















0.1 


0.7 


8.8 


48 

















0.2 


0.3 


1.9 


27 


195 












6 











0.2 


1.3 


8.0 




















0.2 


1.3 


24 


144 














0.2 


0.1 


0.4 


2.7 


74 


544 












7 











0.5 


5.9 






















0.3 


2.1 


53 
















0.2 


0.1 


0.6 


5.0 


152 














8 








0.1 


0.4 


11 






















0.3 


3.2 


72 
















0.2 


0.2 


0.6 


8.5 


216 














10 








0.2 


1.5 





















0.1 


0.4 


6.1 


















0.3 


0.2 


0.8 


12 
















15 








0.3 


3.8 


















0.1 


0.1 


1.2 


17 


















0.2 


0.3 


1.8 


41 
















20 


0.1 


0.1 


0.6 


4.4 


















0.1 


0.3 


2.2 


41 


















0.2 


0.4 


2.9 


73 
















30 


0.1 


0.2 


2.7 


37 


















0.2 


0.5 


5.1 


106 


















0.3 


0.6 


6.8 


170 
















40 


0.3 


0.3 


5.2 


93 


















0.4 


1.1 


10 


242 


















0.6 


1.3 


12 


414 
















50 


0.5 


0.7 


8.2 




















0.6 


1.8 


17 




















0.8 


2.0 


20 
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5.2. Integration of a monomial, using iterated Laurent series. 

In this section, we test the implementation of the method of iterated 
Laurent expansion described in Remark 15 of Section 3.4.'^ 
Table 7 shows the results. 

Table 7. Integration of a random monomial of pre- 
scribed degree using iterated Laurent series. 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 


100 


200 


300 


2 





























0.3 


0.2 






















0.1 


0.1 


0.5 


3.0 


8.1 







0.1 


0.1 


0.1 


0.1 


0.1 


0.2 


0.3 


1.2 


7.1 


30 


3 




















0.1 


0.1 


0.2 

















0.1 


0.4 


0.5 


1.2 


3.4 


34 














0.1 


0.2 


6.2 


3.7 


7.1 


11 


164 






4 














0.1 


0.3 


0.4 


1.0 
















0.1 


0.3 


1.4 


4.8 


16 


39 













0.1 


0.2 


0.7 


12 


35 


77 


176 








5 








0.1 


0.1 


0.1 


1.4 


0.2 












0.1 


0.1 


0.2 


0.7 


4.5 


36 


78 












0.1 


0.1 


0.4 


12 


35 


551 


353 










6 


0.1 





0.1 


0.2 


0.4 
















0.1 


0.2 


0.3 


1.6 


24 
















0.1 


5.0 


0.6 


5.9 


205 














7 


0.1 


0.1 


0.2 


0.4 


















0.1 


0.3 


0.7 


4.1 


















0.2 


5.4 


4.6 


22 
















8 


0.2 


0.2 


0.2 


0.3 


















0.2 


0.3 


1.0 


11 


















0.2 


4.8 


4.1 


111 
















10 


0.3 


0.4 


0.4 




















0.5 


0.5 


2.8 




















7.0 


3.5 


12 


















15 


1.3 


1.4 


3.5 




















1.7 


1.9 


25 




















5.8 


5.1 


73 


















20 


3.8 


4.0 


4.7 




















5.0 


5.3 


123 




















9.4 


8.1 


352 



















30 25 26 
29 29 

41 32 



40 88 90 

98 101 

106 152 

50 248 259 

271 283 

300 429 



®This method is implemented in the Maple procedure integraLvia_iterated, de- 
fined in the file iterated_laurent .mpl. 
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5.3. Integration of a monomial, using Taylor expansion. Here, 
we test the implementation of the algorithm described in Section 3.4. 
This algorithm is based on Corollary 11 
The running times are shown in Table 8. 

Table 8. Integration of a random monomial of pre- 
scribed degree using Taylor expansion. 



n 












Degree 










1 


2 


5 


10 


20 


30 40 


50 


100 


200 


300 


2 


























0.9 


1.9 



















0.1 0.2 


0.3 


2.2 


24 


68 













0.1 


0.1 


0.3 0.6 


0.9 


13 


101 


426 


3 

















0.1 






















0.1 


2.0 


8.1 39 


61 
















0.1 


0.4 


12 


60 277 


512 








4 


































0.1 


1.5 






















0.4 


34 














5 











0.1 






















0.1 


4.9 






















0.6 


48 














6 











0.1 






















0.4 


29 






















1.7 


236 














7 

































0.5 
























1.7 
















8 

































1.1 
























16 
















10 

































3.0 





















0.1 


33 
















15 








0.1 





















0.1 


20 





















0.1 


64 
















20 








1.2 





















0.7 


81 





















28 


205 

















30 

0.7 

0.1 15 

40 

0.1 1.2 

0.1 23 

50 0.1 0.1 

0.1 1.7 

0.2 18 



^This method is implemented in the Maple procedure integraLvia_duality, defined 
in the file duality .mpl. 
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5.4. Integration of a dense homogeneous polynomial over a 
simplex. Following the tests on single monomials, we ran tests on 
random polynomials of varying density. We generated these polyno- 
mials using the Maple function randpoly, requesting a number r of 
monomials and the homogeneous degree M. For each monomial, the 
exponent vector was drawn uniformly from {MgN'^: |M| = M}, and 
the coefficient is drawn uniformly from {!,..., 100}. Due to collisions, 
the generated polynomial actually can have fewer monomials than r.^ 

We only include the results for a family of randomly generated, very 
dense homogeneous polynomials, where we draw r = (^^]~^) random 
monomials. Table 9 shows the number of monomials in the resulting 
polynomials for our random tests. Tables 10, 12, 11 show the test 
results of the three methods. 

We remark that in the case of the method using decompositions 
into powers of linear forms, we note that the same powers of linear 
forms appear in the decomposition formulas (13) for many different 
monomials x^^ , x^^ _ \y"g ^a]^g advantage of this fact by collecting the 
coefficients of powers of linear forms. ^ 

5.5. Integration of a monomial with few effective variables. 

Finally, we tested the performance of the three algorithms on mono- 
mials x^ with a small number D of effective variables. We fix the 
number D. Then, for a given dimension n > D and total degree M, 
we picked 50 combinations of a random simplex A of dimension n and 
a random exponent vector M = (Mi, . . . , Md, 0, . . . , 0) with |M| = M. 
We only include the results for D = 2 in Tables 13, 15, and 14. 



This is implemented in the Maple procedure rctndom_sparse_homogeneous_ 
polynomial_with_degree in the file excunples .mpl. 

^This is implemented in the procedure list_integral_via_waring. 
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Table 9. Number of monomials in a test family of ran- 
dom dense homogeneous polynomials of prescribed de- 
gree 



Degree 





1 


2 


5 


10 


20 


30 


40 


50 


100 


200 


300 


1000 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 




1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 




1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


2 


2 


1 


2 


4 


11 


16 


22 


28 


59 


120 


179 


606 




2 


2.2 


4.2 


7.1 


14 


20 


26 


33 


64 


128 


189 


634 




2 


3 


6 


9 


18 


25 


33 


37 


70 


140 


198 


659 


3 


3 


2 


11 


37 


133 


299 


524 


813 


3209 


12755 








3 


4 


13 


42 


146 


313 


543 


839 


3262 


12836 








3 


6 


17 


46 


160 


325 


560 


867 


3307 


12923 






4 


4 


4 


31 


168 


1085 


3396 


7748 


14719 












4 


6.5 


35 


182 


1118 


3452 


7802 


14807 












4 


8 


40 


196 


1144 


3499 


7887 


14914 










5 


5 


8 


72 


610 


6650 


29181 
















5 


9.7 


81 


633 


6716 


29329 
















5 


12 


88 


661 


6789 


29536 














6 


6 


10 


146 


1847 


33406 


















6 


14 


160 


1896 


33591 


















6 


16 


174 


1934 


33880 
















7 


7 


15 


276 


5015 




















7 


18 


291 


5062 




















7 


21 


304 


5114 


















8 


8 


19 


476 


12213 




















8 


23 


498 


12295 




















8 


29 


521 


12396 


















10 


10 


29 


1236 






















10 


36 


1262 






















10 


41 


1289 




















15 


15 


64 


7274 






















15 


76 


7352 






















15 


85 


7442 




















20 


20 


123 


26742 






















20 


133 


26880 






















20 


146 


26985 




















30 


30 


282 
























30 


295 
























30 


309 























40 40 502 
40 521 

40 542 
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Table 10. Integration of a random dense homogeneous 
polynomial of prescribed degree by decomposition into a 
sum of powers of linear forms 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 100 200 300 


2 











0.1 


1.0 


3.1 


7.8 


16 228 













0.1 


1.2 


3.6 


9.0 


18 270 










0.1 


0.2 


1.4 


4.4 


11 


21 304 


3 








0.1 


2.1 


74 
















0.2 


2.5 


81 










0.1 


0.1 


0.3 


2.9 


88 








4 








0.6 


30 















0.1 


0.8 


32 












0.1 


0.1 


1.0 


35 










5 





0.1 


2.6 

















0.1 


3.0 














0.2 


0.2 


3.4 












6 





0.1 


8.9 














0.1 


0.2 


9.5 














0.2 


0.3 


10 












7 


0.1 


0.2 


26 














0.1 


0.3 


28 














0.2 


0.4 


30 












8 


0.1 


0.4 


81 














0.1 


0.5 


85 














0.2 


0.7 


91 












10 


0.2 


0.9 
















0.3 


1.0 
















0.4 


1.2 














15 


0.6 


4.0 
















0.7 


4.3 
















0.8 


4.8 














20 


1.6 


12 
















1.8 


13 
















1.9 


14 














30 


5.8 


63 
















6.1 


68 
















6.4 


72 














40 


15 


221 















16 232 

16 241 

50 29 
30 

30 
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Table 1 1 . Integration of a random dense homogeneous 
polynomial of prescribed degree using iterated Laurent 
series 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 


100 200 300 


2 














0.2 


0.6 


1.2 


2.1 


24 255 













0.1 


0.3 


0.8 


1.7 


3.1 


29 350 




0.1 


0.1 


0.1 


0.2 


0.4 


1.0 


2.2 


4.0 


37 438 


3 








0.2 


2.1 


22 


146 













0.1 


0.4 


2.7 


32 


184 










0.2 


0.2 


0.6 


3.3 


37 


210 








4 





0.1 


2.4 


37 














0.1 


0.2 


3.0 


40 














0.2 


0.3 


3.6 


44 












5 


0.2 


0.3 


12 


351 














0.2 


0.5 


14 


365 














0.3 


0.7 


16 


380 












6 


0.3 


1.0 


50 
















0.4 


1.2 


55 
















0.5 


1.5 


61 














7 


0.7 


2.0 


161 
















0.8 


2.4 


179 
















0.9 


2.8 


195 














8 


1.2 


3.3 


481 
















1.3 


4.6 


517 
















1.6 


5.8 


560 















10 3.5 13 
3.6 15 

3.8 18 



15 22 117 
23 127 

23 140 



20 91 
94 
98 
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Table 12. Integration of a random dense homogeneous 
polynomial of prescribed degree using Taylor expansion. 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 


100 


2 














0.1 


0.3 


1.2 


1.2 


15 













0.1 


0.2 


0.8 


2.3 


5.6 


98 










0.1 


0.1 


0.3 


1.1 


3.1 


7.4 


129 


3 








0.1 


1.1 


60 


















0.3 


2.4 


132 


















2.0 


3.3 


158 










4 








1.1 


62 




















1.6 


117 

















0.1 


2.6 


143 












r; 



n 
u 


n 1 

U. i 


A A 



















0.1 


8.1 



















0.1 


24 














6 





0.1 


28 



















0.2 


35 
















0.1 


2.3 


41 














7 





0.2 


117 



















0.4 


136 
















0.1 


5.5 


150 














8 





0.4 


366 
















0.1 


0.7 


431 
















0.1 


12 


475 















200 300 



10 0.1 0.8 

0.1 1.3 

0.1 9.8 

15 0.1 4.0 

0.2 5.6 

0.2 7.3 

20 0.3 16 

0.4 18 

2.8 20 

30 0.9 95 

1.1 104 

4.3 109 

40 2.3 370 

2.6 384 

4.4 403 

50 5.2 
5.3 
5.5 



32 V. BALDONI, N. BERLINE, J. A. DE LOERA, M. KOPPE, AND M. VERGNE 



Table 13. Integration of a monomial of prescribed de- 
gree with 2 effective variables by decomposition into a 
sum of powers of linear forms 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 


100 


200 


300 


3 


Q 























1.9 


1.9 


9.8 













0.1 


0.3 


0.7 


1.5 


2.2 


IX 


47 


109 







0.2 


0.1 


0.2 


0.5 


1.1 


2.3 


3.5 


15 


80 


217 


4 























0.5 


0.1 


0.2 


0.3 













0.1 


0.4 


0.9 


1.7 


2.8 


13 


64 


180 







0.3 


0.2 


0.2 


0.7 


1.4 


2.8 


4.3 


20 


103 


270 


5 


























1.8 


0.3 


6.0 










0.1 


0.1 


0.6 


1.3 


2.5 


3.5 


17 


81 


205 







0.4 


0.2 


0.3 


1.0 


2.0 


3.8 


5.9 


26 


128 


335 


g 


























2.1 


7.0 


0.3 










0.1 


0.2 


0.7 


1.8 


3.1 


4.5 


25 


105 


241 




0.2 


0.2 


0.2 


0.4 


1.2 


2.6 


4.9 


7.5 


34 


164 


449 


7 


























2.3 


0.2 


0.2 










0.1 


0.2 


1.0 


1.9 


3.7 


5.9 


24 


128 


279 







0.3 


0.3 


0.5 


1.5 


3.0 


5.9 


8.7 


39 


184 


468 


8 























0.1 


4.2 


0.4 


10.0 










0.1 


0.3 


1.1 


2.3 


4.0 


6.8 


28 


122 


347 







0.5 


0.3 


0.6 


1.6 


3.5 


6.5 


10.0 


45 


206 


529 


10 























1.7 


0.1 


13 












0.1 


0.4 


1.3 


3.3 


5.6 


10 


41 


189 






u.o 


u.o 


n A 


n 7 

U. 1 


2 2 


A Q 




14 


oo 






15 














0.1 


0.1 


0.1 


0.1 


0.2 


0.2 






0.1 


0.1 


0.2 


0.7 


2.4 


6.0 


9.9 


15 


65 


292 






0.3 


0.3 


0.5 


1.2 


4.0 


8.9 


16 


24 


105 


479 




20 





0.1 





0.1 


0.1 


0.1 


0.1 


0.2 


23 








0.1 


0.2 


0.5 


1.3 


4.7 


9.7 


17 


26 


123 








0.3 


0.5 


0.8 


2.1 


7.1 


15 


26 


40 


170 






30 


0.1 


0.2 


0.2 


0.2 


0.2 


0.2 


0.2 


0.2 


0.3 








0.2 


0.4 


0.9 


2.8 


11 


19 


38 


55 


254 








0.4 


0.8 


1.5 


4.7 


15 


33 


58 


92 


369 






40 


0.3 


0.3 


0.3 


0.3 


0.3 


0.4 


0.4 


0.4 










0.4 


0.6 


1.9 


5.4 


20 


41 


66 


101 










0.6 


1.2 


2.9 


8.5 


29 


62 


110 


171 









50 0.5 0.5 0.5 0.6 0.6 0.6 0.7 36 
0.6 1.0 3.4 9.2 31 63 130 201 

0.8 2.1 4.7 14 49 106 185 286 
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Table 14. Integration of a monomial of prescribed 
degree with 2 effective variables using iterated Laurent 
expansion 



Degree 



n 


1 


2 


5 


10 


20 


30 


40 


50 


100 


200 


300 


3 


























0.2 


0.3 


0.9 
















0.1 


0.2 


0.3 


0.3 


1.2 


6.2 


16 













0.1 


0.1 


0.3 


5.6 


2.1 


3.6 


20 


62 


4 














0.1 


0.1 


0.1 


0.1 


0.2 


0.4 


0.9 










0.1 


0.1 


0.1 


0.3 


0.3 


0.5 


2.1 


11 


47 




0.1 


0.1 


0.1 


0.1 


0.2 


4.2 


2.5 


3.0 


5.4 


38 


166 


5 


0.1 





0.1 


0.1 


0.1 


0.1 


0.1 


0.1 


0.3 


0.6 


1.3 




0.1 


0.1 


0.1 


0.1 


0.3 


0.3 


0.5 


0.7 


3.6 


21 


86 




0.1 


0.1 


0.1 


0.2 


6.4 


2.3 


2.8 


3.7 


9.8 


80 


271 


6 


0.1 


0.1 


0.1 


0.1 


0.1 


0.1 


0.2 


0.2 


0.4 


1.0 


1.9 




0.1 


0.1 


0.1 


0.2 


0.3 


0.6 


0.8 


1.1 


7.2 


39 


100 




0.2 


0.1 


0.2 


3.4 


2.5 


2.6 


3.8 


4.7 


30 


269 


550 


7 


0.1 


0.1 


0.1 


0.1 


0.2 


0.2 


0.2 


0.3 


0.6 


1.2 






0.1 


0.1 


0.2 


0.3 


0.6 


0.8 


1.2 


2.0 


6.8 


94 






0.2 


0.2 


3.9 


2.5 


2.8 


3.4 


4.2 


6.7 


26 


260 




8 


0.2 


0.2 


0.2 


0.2 


0.2 


0.2 


0.3 


0.3 


0.7 


2.0 






0.2 


0.3 


0.3 


0.4 


0.7 


1.2 


1.9 


2.9 


16 


99 






0.3 


4.5 


2.3 


2.6 


3.1 


3.6 


5.8 


12 


42 


322 




10 


0.3 


0.3 


0.3 


0.4 


0.4 


0.5 


0.6 


0.6 


1.1 








0.5 


0.4 


0.5 


0.7 


1.2 


2.4 


3.8 


5.9 


23 








4.2 


2.3 


2.5 


3.0 


3.9 


6.5 


11 


17 


77 






15 


1.3 


1.3 


1.4 


1.4 


1.6 


1.7 


1.9 


2.4 










1.4 


1.5 


1.7 


2.1 


4.0 


7.3 


11 


20 










1.6 


1.9 


2.0 


3.4 


9.1 


21 


31 


52 








20 


4.0 


4.0 


4.1 


4.2 


















4.2 


4.4 


4.8 


6.1 


















4.4 


5.6 


5.7 


9.3 
















30 


22 


24 


24 


25 


26 
















25 


25 


27 


32 


47 
















26 


28 


31 


45 


96 














40 


71 


88 


90 




















85 


92 


95 




















92 


102 


103 


















50 


196 


228 


212 


217 


















231 


239 


244 


261 


















252 


249 


277 


347 
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Table 15. Integration of a monomial of prescribed de- 
gree with 2 effective variables using Taylor expansion 



Degree 







2 


5 


10 


20 


30 


40 


50 


100 


900 '^00 


3 





























1.9 
















0.1 


0.2 


0.4 


0.8 


5.8 


62 










0.1 


0.1 


0.2 


0.6 


1.3 


2.6 


27 


359 


4 













































0.2 


0.4 


0.7 


1.8 


11 












0.1 


0.1 


0.4 


1.2 


3.0 


8.2 


58 




5 










































0.1 


0.2 


0.6 


2.0 


2.9 


19 












0.1 


0.2 


0.7 


2.2 


6.7 


15 


110 




6 










































0.1 


0.7 


0.9 


2.0 


4.4 


31 












0.1 


0.3 


16 


8.4 


8.8 


72 


124 




7 




















0.2 





















0.1 


0.5 


1.6 


4.3 


5.9 


52 









0.1 


0.1 


0.3 


1.8 


8.4 


21 


32 


246 




8 























0.3 


0.1 












0.1 


0.2 


0.9 


2.7 


3.6 


9.3 


65 






0.1 


0.1 


0.2 


0.4 


2.8 


36 


19 


42 


547 




10 






































0.1 


0.2 


1.1 


3.5 


8.4 


20 








0.1 


0.1 


0.2 


0.6 


7.0 


22 


40 


92 






20 























2.0 











0.1 


0.4 


2.8 


12 


40 


74 


121 








0.1 


0.2 


1.1 


9.3 


43 


156 


309 


462 






30 

















0.1 















0.2 


1.4 


6.9 


62 


90 












0.1 


0.5 


5.6 


29 


205 


583 










40 








0.1 


0.1 


0.1 














0.1 


0.4 


3.3 


10 


71 














0.1 


0.9 


32 


40 


329 












50 


0.1 


0.1 


0.1 


0.1 


0.1 














0.1 


0.4 


5.5 


19 


126 














0.1 


1.5 


39 


65 


525 
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5.6. Discussion. In our implementation of the three methods and our 
experiments for the case of random monomials, we observe that the 
method of iterated Laurent expansion is faster than the two other 
methods if the dimension n is very small (up io n = 5). Starting 
from dimension n = 6, the method using decompositions into powers 
of linear forms is faster than the other two methods. The method us- 
ing Taylor expansion is always inferior to the better of the two other 
methods, for any combination of degree and dimension. 

In the experiments with random dense polynomials, in our implemen- 
tation we did not see significant savings from collecting the coefficients 
of the same powers of linear forms. As a consequence, the ranking of 
the three methods is the same as it is in the case of random monomials. 

The experiments with random monomials with few effective variables 
show that all three methods benefit from using few effective variables. 
The greatest effect is on the method using decompositions into pow- 
ers of linear forms, where, for example, the restriction to 2 effective 
variables allows to handle combinations of high degree M = 200 and 
high dimension n = 15. However, for low dimensions (n < 5), the 
method of iterated Laurent expansion still wins. Also here the method 
using Taylor expansion is always inferior to the better of the two other 
methods. This discussion shows the power of Brion's formula. 
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6. Conclusions 

We discussed various algorithms for the exact integration of polyno- 
mials over simplicial regions. Beside their theoretical efficiency, the sim- 
ple rough experiments we performed clearly demonstrated that these 
methods are robust enough to attack rather difficult problems. Our 
investigations opened several doors for further development, which we 
will present in a forthcoming paper. 

First, we have some theoretical issues expanding on our results. As 
in the case of volumes and the computation of centroids, it is likely that 
our hardness result. Theorem 1, can be extended into an inapproxima- 
bility result as those obtained in [28]. Another goal is to study other 
families of polytopes for which exact integration can be done efficiently. 
Furthermore, we will present a natural extension of the computation of 
integrals, the efficient computation of the highest degree coefficients of 
a weighted Ehrhart quasipolynomial of a simplex. Besides the methods 
of the present article, these last computations are based on the results 
of [7] and [6]. 

Second, our intention has been all along to develop algorithms with 
a good chance of becoming practical and that allow for clear imple- 
mentation. Thus we have also some practical improvements to discuss. 

Finally, in order to develop practical integration software, it appears 
that our methods should be coupled with fast techniques for decom- 
posing domains into polyhedral regions (e.g. triangulations) . 
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